# load the data↔
I observed that the price of options on the day that they expire is odd. If one bought or sold options on the day of their expiry one could frequently make profits:
# find arbitrage opportunities (buy an option and exercise it, or sell an option and fulfill it)↔
# plot histogram of Dollar profits↔
I conclude that I will never want to use the price of an option on the day of its expiry. Hence, in (b) I choose a minimal duration to be a few days longer than the desired holding period.
On every day for which there is options data available I sell one ATM call and and one ATM put and receive the best_bid price. I then look up the best_offer price of the respective options via optionid on the first day after the holding period ends.
I use two different schemes to select the options that are sold, and then bought back in a month:
Note: The way the return is calculated is very important!
At first, I calculate returns as if the investor did not have to adhere to margin requirements. That is, there is no necessary deposit at a margin account. Hence, the formula for the return is:
$$R^\text{naive}_{t,t+1} = \frac{P^{\text{call}, i}_t + P^{\text{put}, j}_t}{P^{\text{call}, i}_{t+1} + P^{\text{put}, j}_{t+1}} - 1$$
Then, at the end of (b), (c), and (d) I report returns for a strategy where the margin account needs to have a balance of $M_t = \Big(P^{\text{call}, i}_t + P^{\text{put}, j}_t\Big)\Big(\frac{1}{1-0.8} - 1\Big)$. That is, the the balance on the margin account plus the proceeds from selling the options needs to be enough to cover five times the the current option price. This is motivated by the maximal observed loss of the strategy of $-0.7946$. The return is then:
$$R^\text{margin}_{t,t+1} = \frac{P^{\text{call}, i}_t + P^{\text{put}, j}_t - P^{\text{call}, i}_{t+1} - P^{\text{put}, j}_{t+1}}{M_t} + R^\text{risk free}_t$$
function ATM_focused_scheme(options, trading_dates; last_day = Date(2016, 3, 31), min_days = 33)↔
end
function TTM_focused_scheme(options, trading_dates; last_day = Date(2016, 3, 31), min_days = 33)↔
end
I evaluate both schemes:
returns_atm = ATM_focused_scheme(options, trading_dates);
returns_ttm = TTM_focused_scheme(options, trading_dates);
There are not many errors keeping in mind that I evaluate the strategy for every day.
The ATM-focused scheme is more concentrated around the zero percent deviation from ATM, while the TTM-focused scheme has shorter time-to-maturity. The holding periods are equal. And the two scheme regularly select different options:
# plots to assess how well the strategy is fulfilled↔
I execute the trading strategy on the end of every month. The two schemes perform very similarly:
# plot the performance↔
The returns of the ATM-focused scheme with a four-fold margin requirement are much more realistic:
# plot the performance↔
returns_atm1y = ATM_focused_scheme(options, trading_dates, last_day = Date(2015, 3, 31), min_days = 368)
returns_ttm1y = TTM_focused_scheme(options, trading_dates, last_day = Date(2015, 3, 31), min_days = 368);
There are some more instances where there is no price at the desired buy-back date and quite a few times there are not enough contracts traded on a sell day.
Again, the ATM-focused scheme yields lower deviation from ATM, while the TTM-focused scheme has no time-to-maturity outliers but has moneyness outliers. The strategies agree quite often, but with episodes of significant deviation:
# plots to assess how well the strategy is fulfilled↔
If I chose to execute the strategy on the last trading date of every year, the performance would be:
# plot the performance↔
Let me compare annual returns, that is I report the returns from the monthly strategy started on any given day and repeated 11 times in the following year.
# plot the histogram of annualized returns↔
This histogram suggests that the short-term volatility selling strategy is better than the long-term strategy.
# plot the performance↔
There is no money to be made selling long-term volatility!
# plot the histogram of annualized returns↔
# load risk free rate↔
I use the ATM-focused excess returns and calculate point estimates. The monthly returns are aggregate to yearly returns:
returns1m = ([prod([returns_atm[searchsortedlast(returns_atm[:date], d + Dates.Month(m)), :return] for m in 0:11] .+ 1) - 1 for d in returns_atm[:date][1:4843]]+1) ./ ([interest[searchsortedlast(interest[:date], d), :RF] for d in returns_atm[1:4843, :date]]/100 + 1) .- 1
returns1y = ((returns_atm1y[:return]+1) ./ ([interest[searchsortedlast(interest[:date], d), :RF] for d in returns_atm1y[:date]]/100 + 1) .- 1)[1:4843];
# point estimates with standard functions↔
But I need an approximate distribution, so I do joint GMM (i.e. both the monthly and yearly strategy in one estimation, this turns out not to be relevant) with the following moment function, where $\theta_1 = \mu_{1m}$, $\theta_2 = \sigma_{1m}$, $\theta_1 = \mu_{1y}$ and $\theta_2 = \sigma_{1y}$:
#operates on a single observation, w[1] = return_1m, w[2] = return_1y
function g(w, θ)
[w[1] - θ[1], (w[1] - θ[1])^2 - θ[2]^2, w[2] - θ[3], (w[2] - θ[3])^2 - θ[4]^2]
end
# my standard GMM setup, with Newey-West HAC variance estimation↔
# standard efficient GMM procedure↔
result = eff_GMM(hcat(returns1m, returns1y), Jn = 30)
This yields: $$\left[\begin{array}{c} \mu_{1m}\\ \sigma_{1m}\\ \mu_{1y}\\ \sigma_{1y} \end{array}\right]\sim\mathcal{N}\left(\left[\begin{array}{c} 12.2006\\ 37.8966\\ 0.19194\\ 0.69552 \end{array}\right],\left[\begin{array}{cccc} 2.4026 & 9.18919 & 0.029511 & 0.041627\\ 9.18919 & 44.6641 & 0.11467 & 0.210015\\ 0.029511 & 0.11467 & 0.002228 & 0.002359\\ 0.041627 & 0.21005 & 0.002359 & 0.005462 \end{array}\right]\right)$$
Applying the Delta-Method with $f(\theta) = \big[\frac{\theta_1}{\theta_2}, \, \frac{\theta_3}{\theta_4}\big]$:
f(θ) = [θ[1]/θ[2], θ[3]/θ[4]]
SR_result = [f(result[1]), ForwardDiff.jacobian(f, result[1]) * result[2] * ForwardDiff.jacobian(f, result[1])']
That is:
$$\left[\begin{array}{c} SR_{1m}\\ SR_{1y} \end{array}\right]\sim\mathcal{N}\left(\left[\begin{array}{c} 0.32194\\ 0.27597 \end{array}\right],\left[\begin{array}{cc} 0.0007764 & -0.000008\\ -0.000008 & 0.002775 \end{array}\right]\right)$$
Hence, the difference $f(\theta)=\theta_1-\theta_2$:
f(θ) = [θ[1] - θ[2]]
Diff_result = [f(SR_result[1]), ForwardDiff.jacobian(f, SR_result[1]) * SR_result[2] * ForwardDiff.jacobian(f, SR_result[1])']
That is:
$$\Delta SR = SR_{1m}-SR_{1y}\sim\mathcal{N}\left(0.04597,\,0.00357\right)$$
The $\mathcal{H}_0$ is that $\Delta SR \leq 0$, the $\mathcal{H}_1$ is that $\Delta SR \gt 0$. The test statistic is $\frac{\mu_{\Delta SR}}{\sigma_{\Delta SR}}$, and the corresponding p-value can be calculated by see how much mass there is below zero:
using Distributions
cdf(Normal(Diff_result[1][1], Diff_result[2][1]), 0.)
A clear rejection of the $\mathcal{H}_0$!
# point estimates with standard functions↔
# running GMM and the Delta method again, evalutating the test statistic↔
The p-value is again significant! The $\mathcal{H}_0$ is rejected!
# load the data↔
My convention is that the date in the row marks the end-of-period, the return and rf has been achieved inside of this period, but the vrp only became available at the end. Multi-period returns start in the period for with Date is the end-of-month date.
Hence in predictive regressions I need to regress next period's return on this period's vrp.
function OLS(data, sym)
df = dropmissing(data[[sym, :vrp]])
Y = convert.(Float64, df[2:end, sym])
X = hcat(ones(length(Y)), df[1:end-1, :vrp])
β = round.((X'X) \ (X'Y), 5)
ϵ = Y - X * β
σsq = (ϵ' * ϵ) / length(ϵ)
aVarβ = (X' * X) ./ length(ϵ) .* σsq
DataFrame(LHS = string(sym), End_Date = data[end, :Date], α = β[1], α_sd = sqrt.(diag(aVarβ) ./ length(ϵ))[1], β = β[2], β_sd = sqrt.(diag(aVarβ) ./ length(ϵ))[2])
end
#tabular(vcat(vec([OLS(data[ran, :], Symbol("ex_return1$p")) for p in ["m", "q", "y"], ran in [1:216, 1:size(data, 1)]])...), rounding=4);
\begin{array}{cccccc} \hline LHS & until & \beta_0 & \sigma(\beta_0) & \beta_1 & \sigma(\beta_1)\\ \hline \text{1m} & 2007 & -0.0002 & 0.0026 & 0.0003 & 0.0637\\ \text{1q} & 2007 & -0.0064 & 0.0044 & 0.0013 & 0.106\\ \text{1y} & 2007 & 0.0323 & 0.0102 & 0.002 & 0.2466\\ \text{1m} & 2017 & -0.0007 & 0.0022 & 0.0004 & 0.0565\\ \text{1q} & 2017 & 0.0017 & 0.0038 & 0.0011 & 0.0984\\ \text{1y} & 2017 & 0.0645 & 0.0088 & 0.0012 & 0.2318 \end{array}
After the financial crisis the predictive power of the variance risk premium is getting weaker using it for longer horizons.
Mind: the reported standard errors are not adjusted for autocorrelation!
function GLS(data, sym)
df = dropmissing(data[[sym, :vrp, :vix]])
Y = convert.(Float64, df[2:end, sym])
X = hcat(ones(length(Y)), df[1:end-1, :vrp])
β = round.((X' * diagm(1 ./ df[1:end-1, :vix]) * X) \ (X' * diagm(1 ./ df[1:end-1, :vix]) * Y), 5)
ϵ = Y - X * β
σsq = (ϵ' * ϵ) / length(ϵ)
aVarβ = (X' * X) ./ length(ϵ) .* σsq
DataFrame(LHS = string(sym), End_Date = data[end, :Date], α = β[1], α_sd = sqrt.(diag(aVarβ) ./ length(ϵ))[1], β = β[2], β_sd = sqrt.(diag(aVarβ) ./ length(ϵ))[2])
end
#tabular(vcat(vec([GLS(data[ran, :], Symbol("ex_return1$p")) for p in ["m", "q", "y"], ran in [1:216, 1:size(data, 1)]])...), rounding=4);
\begin{array}{cccccc} \hline LHS & until & \beta_0 & \sigma(\beta_0) & \beta_1 & \sigma(\beta_1)\\ \hline \text{1m} & 2007 & 0.0018 & 0.0026 & 0.0002 & 0.0638\\ \text{1q} & 2007 & 0.0008 & 0.0044 & 0.001 & 0.1063\\ \text{1y} & 2007 & 0.0478 & 0.0102 & 0.0016 & 0.2472\\ \text{1m} & 2017 & 0.001 & 0.0022 & 0.0004 & 0.0565\\ \text{1q} & 2017 & 0.0053 & 0.0038 & 0.001 & 0.0985\\ \text{1y} & 2017 & 0.0683 & 0.0088 & 0.0012 & 0.2319 \end{array}
Bollerslev et al (2009) define $VRP_t = IV_t - RV_t$ where $RV_t$ is the realized volatility in the period that ended with $t$. $IV_t$ is the implied volatility.
Hence VIX $\approx IV_t$, because the VIX is the implied volatility from one-month options. GLS down-weights exactly those observations where $IV_t$ is high. This decreases the "freak" losses in the financial crisis, when the VIX stayed high for some time.
So the GLS estimation predicts positive returns whenever the VIX is high and the realized volatility is low, but down-weighs observations where VIX is high. In sum, observations with low realized volatility are over-weighted.
This makes sense if I think of the variance risk premium as being paid if the volatility stays low after being low already.